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Abstract 

The equations of motion of compact binary systems have been derived in the post-Newtonian (PN) approximation 
of general relativity. The current level of accuracy is 3.5PN order. The conservative part of the equations of 
motion (neglecting the radiation reaction damping terms) is deducible from a generalized Lagrangian in harmonic 
coordinates, or equivalently from an ordinary Hamiltonian in ADM coordinates. As an application we investigate 
the problem of the dynamical stability of circular binary orbits against gravitational perturbations up to the 3PN 
order. We find that there is no innermost stable circular orbit or ISCO at the 3PN order for equal masses. 

Resume 

Les equations du mouvement d'un systeme binaire d'objets compacts out etees calculees dans I'approximation 
post-newtonienne (PN) de la relativite generale. Le niveau d'approximation atteint I'ordre 3.5PN. La partie 
conservative des equations du mouvement (obtenue en negligeant les termes de freinage de rayonnement) se deduit 
d'un lagrangien generalise en coordonnees harmoniques, et, de fagon equivalente, d'un hamiltonien ordinaire 
en coordonnees ADM. Comme application nous etudions le probleme de la stabilite dynamique, vis-a-vis des 
perturbations gravitationnelles, des orbites binaires circulaires a I'ordre 3PN. Nous trouvons qu'il n'y a pas de 
derniere orbite stable circulaire on ISCO a I'ordre 3PN pour des masses egales. 
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1. Introduction 

The problem of the dynamics of two compact bodies is part of a program aimed at unravelling the information 
contained in the gravitational-wave signals emitted by inspiralling and coalescing compact binaries — systems of 
neutron stars or black holes driven into coalescence by emission of gravitational radiation. The treatment of the 
problem is Post-Newtonian, i.e. based on formal expansions, when the speed of light c tends to infinity, of general 
relativity theory around Newton's theory. The early, classic works of Lorentz & Droste [1], Eddington & Clark [2], 
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Einstein, Infeld & Hoffmann [3] , Fock [4] , Papapetrou [5] and others led to a good understanding of the equations 
of motion of A'^ bodies at the first post-Newtonian approximation (IPN, corresponding to order ~ In the 

1970's, a series of works [6,7,8] led to a nearly complete control of the problem of motion at the second post- 
Newtonian approximation (2PN ^ 1/c''). Then, starting in the early 80's, motivated by the observation of secular 
orbital effects in the Hulse- Taylor binary pulsar PSR1913+16 [9,10], several groups solved the two-body problem 
at the 2.5PN level (while completing on the way the derivation of the 2PN dynamics) [11,12,13,14,15,16,17,18]. 
The 2.5PN term constitutes the first contribution of gravitational radiation reaction in the equations of motion 
(such term is the analogue of the Abraham-Lorentz reaction force in electromagnctism) , and is directly responsible 
for the decrease of the binary pulsar orbital period by emission of gravitational radiation. 

In the late 90's, motivated by the aim of deriving high-accuracy templates for the data analysis of the in- 
ternational network of interferometric gravitational-wave detectors LIGO/VIRGO, two groups embarked on the 
derivation of the equations of motion at the third post-Newtonian (3PN) level. One group used the Arnowitt- 
Deser-Misner (ADM) Hamiltonian formalism of general relativity [19,20,21,22] and worked in a corresponding 
ADM-type coordinate system. Another group used a direct post-Newtonian iteration of the equations of motion 
in harmonic coordinates [23,24,25,26,27,28]. The end results of these two approaches have been proved to be 
physically equivalent [22,27]. However, both approaches, even after exploiting all symmetries and pushing their 
Hadamard-regularization-based methods to the maximum, left undetermined one dimensionless parameter at the 
3PN order. The appearance of this unknown parameter was related with the choice of the regularization method 
used to cure the self-field divergencies of point particles. Both lines of works regularized the self-field divergencies 
by some version of the Hadamaxd regularization method. Finally, the completion of the equations of motion at 
the 3PN order was made possible thanks to the powerful dimensional self-field regularization, which could fix up 
uniquely the value of the ambiguity parameter [29,30]. This result is also in complete agreement with the recent 
finding of [31,32], who derived the 3PN equations of motion in harmonic gauge using a "surface-integral" approach 
without use of self-field regularization. Finally, the 3.5PN terms, which constitute a IPN relative modification of 
the radiation reaction force (and are relatively easier to derive), have been added in Refs. [33,34,35,36,37,38]. 

In Section 2 of the present paper we give the final result for the equations of motion of compact binaries at the 
3.5PN order. The equations are presented in ready to use quasi-Newtonian form, in the reference frame associated 
with the center of mass position. In Section 3 we discuss the Lagrangian and Hamiltonian formulations of the 
conservative part of the equations, obtained by neglecting the radiation reaction terms occuring at the 2.5PN 
and 3.5PN orders. Finally, in Section 4, we investigate, following Ref. [28], the question of the stability of circular 
orbits, against linear gravitational perturbations, up to the 3PN order. 



2. Equations of motion in the center-of-mass frame 

In the present paper we employ the so-called haxmonic-coordinates approach [23,24,27,28] which derived the 
3PN binary's equations of motion in harmonic coordinates. From these equations, obtained at first in a general 
frame, we translate the origin of coordinates to the binary's center-of-mass by imposing that the binary's center- 
of-mass vector position is d = 0. The center-of-mass vector is nothing but the conserved integral of the motion 
that is associated, via the Noether theorem, with the boost symmetry of the Lagrangian from which the 3PN 
equations of motion are derived [27,28]. The condition Gj = results in the 3PN-accurate relationship between 
the individual positions of the particles yl and j/2 in the center-of-mass frame, and the relative position a;' = j/i — j/2 
and velocity = v{ — V2 = dx''/dt. The center-of-mass equations of motion are then obtained by replacing the 
individual positions and velocities by their center-of-mass expressions, given in terms of and velocity v^, applying 
as usual the order-reduction of all accelerations where necessary. Order reduction means that any acceleration 
(or derivative of acceleration) in a sub-dominant post-Newtonian term is to be replaced by its explicit expression 
given as a function of the positions and velocities as deduced from the lower-order equations of motion themselves. 

We shall denote the orbital separation by r = |x|, and pose n — :K/r and r = n • v. The mass parameters are 
the total mass m = mi + m2, the mass difference 5m = mi — m2, the reduced mass /i = mimg/m, and the very 
useful symmetric mass ratio v = /x/m, which is such that < < 1/4, with i/ = 1/4 in the case of equal masses, 
and 1/ — » in the "test-mass" limit for one of the bodies. We write the relative acceleration in the center-of-mass 
frame in the form (we pose G = 1) 




(1) 



2 



where the first term represents the famous Newtonian approximation, and where the post-Newtonian remainder 
term 0{c~^) indicates the level of accuracy of the expression which is here 3.5PN order. We find [28,38] that the 
coefficients A and B are 
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The 3.5PN equations of motion play a crucial role when deriving the high-accuracy templates which will be used 
for analysing (hopefully in a near future) the gravitational wave signals from compact binary inspiral in the data 
analysis of the LIGO and VIRGO detectors. 

At the 3PN order we find some logarithmic terms, depending on some arbitrary constant r^. The presence of 
these logarithms reflects in fact the use of a specific harmonic coordinate system. It is indeed known that the 
logarithms at the 3PN order in Eqs. (2), together with the constant therein, can be removed by applying a 
gauge transformation. This shows that there is no physics associated with them, and that these logarithms and 
the constant Tq will never appear in any physical result derived from these equations, because the physical results 
must be gauge invariant (ro is sometimes referred to as a "gauge constant"). The gauge transformation at 3PN 
order whose effect is to remove the logarithms is given in [24] . Notice that after applying this gauge transformation 
we are still within the class of harmonic coordinates. The resulting modification of the equations of motion affects 
only the coefficients of the 3PN order in Eqs. (2); let us denote them by .43pn and S3PN. The new values of 
these coefficients, say ^spn and B'^p^, obtained after removal of the logarithms by the latter harmonic gauge 
transformation, are then [39] 

„ 1 r 35r^V nbi^'u^ 175 r'V'^ l^f^uv^ 

^3PN = ^ { + T7. T7. + 



c6 1^ 16 16 16 2 

_ 135 r" 255 v"^ _\hr^v 237 i/^ 

4 ' 8 2 ^ 8 

45r^i/=S'' Wvv^' 49u^v*' ,.,36 
-^^ + — — + lSr^''v^ 

, m ( .4 69r*i/2 .4 3 .2 2 , .2 2 2 

H 79 r v 30 r v — 121r vv -|-16r v v 

r \ 2 

, ori -2 3 2 , 75z/i;* 2 4 irv 3 4\ 

+ 20r z/ V H : \-8i' v — lOi^ ?; 

4 ) 

(.2 22717 r'^v llr'^u'^ ^.2 3 615r''uw'^ 
H 5- H 1 ^ ' ^' ^ + ■ 



r2 V 168 8 64 

3 2 123 utt'^ 
+ v V — 



840 64 
[ 13991/ 71i/2 411/77^ 
^ — \-^^-—2 2- + ^6- 

w 1/ 45r^^_.5 2^15r^^,„.3 2 

lllf^v'^v'^ 10 -3 3 2 GSrVV^ 2 4 , fi. 34 

12r u V \-19rv v +bru v 

4 8 

, m f 329r^v , 59r'i/' , .3 3 ... 2 . 2 2 . 3 2 

H — 1 |-18r v —Ibrvv —27rv v —lOri/ v 

r \ 6 2 

/ . 5S49rp . 2 „• 3 123rv-K^\] 
+ - (^-4.-^^+25..2 + 8..3 (3) 

These gauge-transformed coefficients are useful because they do not yield the usual complications associated with 
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logarithms. However, they must bo handled with care in applications such as [39], since one must ensure that 
all other quantities in the problem (energy, angular momentum, gravitational- wave fluxes etc.) are defined in the 
same specific harmonic gauge avoiding logarithms. In the following we shall no longer use the coordinate system 

leading to Eqs. (3). Notably the expression wc shall derive below for the Lagrangian will bo valid in the "standard" 
harmonic coordinate system in which the equations of motion are given by (1) with (2). 



3. Lagrangian and Hamiltonian formulations 



The Lagrangian for the relative ccnter-of-mass motion is obtained from the 3PN ccntcr-of-mass equations of 
motion (l)-(2) in which one ignores the radiation-reaction terms at the 2.5PN and 3.5PN orders. We are indeed 
interested in the conservative part of the equations of motion, excluding the terms associated with gravitational 
radiation; only the conservative part is deducible from a Lagrangian. It is known that the Lagrangian in harmonic 
coordinates will necessarily be a generalized one (from the 2PN order), i.e. depending not only on the positions 
and velocities of the particles, but also on their accelerations [13,15]. It is also known [15,27] that one can always 
restrict ourselves to a Lagrangian that is linear in the accelerations. 

The conservative part of the center-of-mass equations of motion (l)-(2) then take the form (after systematic 
order-reduction of the accelerations) of the generalized Lagrange equations 

where L[x^,v^,a^] denotes the generalized center-of-mass Lagrangian — which is linear in the accelerations ap- 
pearing at 2PN and 3PN orders. We recall that there is a large freedom for choosing a Lagrangian because we 
can always add to it the total time derivative of an arbitrary function. As a matter of convenience, we choose 
below a particular ccnter-of-mass Lagrangian that is "close" (in the sense that many coefficients arc identical) 
to some "fictitious" Lagrangian obtained from the general-frame one given in Ref. [27] by the mere Newtonian 
center-of-mass replacements yl Ih — > We point out that such a fictitious Lagrangian is not 

the correct Lagrangian for describing the center-of-mass relative motion. Indeed, the actual relations connecting 
the center-of-mass variables y\ and j/2 to the relative position and velocity u\ involve many post-Newtonian 
corrections, so the actual center-of-mass Lagrangian must contain some extra terms in addition to those of the 
latter fictitious one. However, we find that these extra terms arise only from the 2PN order. Our result (when 
divided by the reduced mass n = mv) is then 
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Witness the acceleration terms present at the 2PN and 3PN orders; our notation is a„ = a ■ n and = a • v for 
the scalar products between a* = dv^ /dt and the direction and velocity n*. We recall here that it is in general 
forbidden to order-reduce the accelerations in a Lagramgian. 

We next consider the problem of the Hamiltonian associated with (the conservative part of) the equations of mo- 
tion (l)-(2). This problem is not straightforward in harmonic coordinates because of the presence of accelerations 
at 2PN and 3PN orders in the Lagrangian (5). To proceed, the best is again to change coordinates, and transform 
the harmonic coordinate system into a new system which avoids the appearance of the accelerations terms, and 
also, as it will turn out, of the logarithms at the 3PN order. This new coordinate change will thus contain a piece 
which is identical to the one we used to remove the logarithms of the equations of motion in Section 2. However, 
this new coordinate system is not harmonic; it was introduced long ago by Arnowitt, Deser & Misner in their study 
of the Hamiltonian formulation of general relativity and is called the ADM coordinate system. In Rcf. [27] the 
"contact" transformation between the particles' variables in harmonic coordinates and those in ADM coordinates 
was determined. By contact transformation we mean the relation between the particles' trajectories y\{t) in one 
coordinate system, and the corresponding trajectories YX{t) in another, say 6yA{t) = YX{t) — j/aC^)- Notice that 
in the contact transformation, the time t is to be viewed as a "dummy" variable. The contact transformation 
is not a coordinate transformation between the spatial vectors in both coordinates, but takes also into account 
the fact that the time coordinate changes as well; i.e. Sy\ = CiVA) — ^°(2/a)i'a/c, where ^''(j/a) denotes the 
four-dimensional change between the coordinates, when evaluated at the position t/a = [t, y\{t)\. 

There is a unique contact transformation such that the 3PN harmonic-coordinates Lagrangian (5) is changed 
into another Lagrangian whose Legendre transform coincides with the 3PN ADM-coordinates Hamiltonian derived 
in [19] (see [27,28] for details). In a first stage this yields the expression for the ADM-coordinates Lagrangian, in 
which we use names appropriate to the ADM variables X"^ = -\- 5x^ , which means the separation distance R, the 
relative square velocity V^, and the radial velocity = N • V. This is an ordinary Lagrangian, depending only on 
the positions and velocities and without accelerations, L^^^[X^, V^], and as we said which is free of logarithms 
at the 3PN order. Its explicit expression is 
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Next wo apply the ordinary Lcgcndrc transform to obtain the corresponding Hamiltonian, H^^^[X^ , P''], which 
is a function of X' and the conjugate momentum = dL'^^^/dV''. We find 
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We denote P^ = and Pr — 'N ■ P. The previous result is in perfect agreement with the center-of-mass 
Hamiltonian derived in Ref. [19]. 



4. Dynamical stability of circular orbits 

As an application let us investigate the problem of the stability, against dynamical gravitational perturbations, of 
circular orbits at the 3PN order. We want in particular to discuss the existence (or non-existence) of an innermost 
stable circular orbit (ISCO) at various post-Newtonian orders, which would constitute the analogue for two black 
holes with comparable masses of the famous orbit risco = 6M/c^ in the Schwarzschild metric. Notice that we are 
concerned here with the stability of the orbit with respect to purely gravitational perturbations appropriate to the 
motion of black holes; however it is known that for neutron stars instead of black holes the ISCO is determined 
by the hydrodynamical instability rather than by the effect of general relativity. 

We propose to use two different methods for this problem, one based on a perturbation at the level of the 
equations of motion (l)-(2) in harmonic coordinates, the other one consisting of perturbing the Hamiltonian 
equations in ADM coordinates for the Hamiltonian (7). We shall find a criterion for the stability of orbits and 
shall present it in the form of an invariant expression (which is the same in different coordinate systems) . We shall 
check that our two methods agree on the result. 

We deal first with the perturbation of the equations of motion, following the approach proposed in Section 111. A 
of Ref. [40]. We introduce polar coordinates (r, <^) in the orbital plane and pose u = r and ui = if. Then Eq. (1) 
yields the system of equations 
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where A and B are given by Eqs. (2) as functions of r, u and lo (through — + r'^uj^). In the case of 
which is circular, apart from the adiabatic inspiral due to the 2.5PN and 3.5PN radiation-reaction effects, 
r = u = ui = hence m = 0. Eq. (8b) gives thereby the angular velocity ujo of the circular orbit as 
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Solving iteratively this relation at the 3PN order using the equations of motion (l)-(2) we obtain wo as a 
of the circular-orbit radius ro in harmonic coordinates (the result agrees with the one of Refs. [23,24]), 
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The circular-orbit radius ro should not be confused with the constant r'o entering the logarithm at the 3PN order 
and which is issued from Eqs. (2). 

Now we investigate the equations of linear perturbations around the circular orbit defined by the constants ro, 
Mo = 0, or, rather, if we were to include the radiation-reaction damping, mo = 0{c~^), and ojo- We pose 

r = ro + 5r , (Ha) 
u = 6u , (lib) 
u! = ujo + 5u! , (11c) 
where Sr, 5u and 5oj denote some perturbations of the circular orbit. Then a system of linear equations follows as 



5r^5u, (12a) 

5u = aa 5r + Po 5ijj , (12b) 

5uj — 70 Su , (12c) 

where the coefficients, which solely depend on the unperturbed circular orbit, read [40] 
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In obtaining Eqs. (13) we use the fact that .4 is a function of the square through = -t-r^w^, so that dA/du 

is proportional to u and thus vanishes in the unperturbed configuration (because u = 5u). On the other hand, 
since the radiation reaction is neglected, B also is proportional to u [see Eqs. (2)], so only dB/du can contribute 
at the zeroth perturbative order. Now by examining the fate of perturbations that are proportional to some e"'*, 
we arrive at the condition for the frequency a of the perturbation to be real, and hence for stable circular orbits 
to exist, as being [40] 

C„ = -ao - /3o 70 > . (14) 
Substituting into this A and B at the 3PN order we then arrive at the orbital-stability criterion 
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where we recall that ro is the radius of the orbit in harmonic coordinates. 

Our second method is to use the Hamiltonian equations based on the 3PN Harriiltonian in ADM coordinates 
given by Eq. (7). We introduce the polar coordinates (-R, *t) in the orbital plane — we assume that the orbital 
plane is equatorial, given by © = f in the spherical coordinate system {R, ©, ^) — and make the substitution 

P' = Pn+^, (16) 

into the Hamiltonian. This yields a "reduced" Hamiltonian that is a function of R, Pr and P*, namely Ti. = 
?^[P, Pr,P*], and describes the motion in polar coordinates in the orbital plane (henceforth we denote H = 
H °'^//Lt). The Hamiltonian equations then read 

(17a) 
(17b) 
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Evidently the constant is nothing but the conserved angulax-momentum integral. For circular orbits we have 
R = Ro (a constant) and Pr = 0, so 

|^[i?o,0,P^] =0, (18) 
which gives the angular momentum of the circular orbit as a function of Ro, and 

which yields the angular frequency of the circular orbit cuo — the same as in Eq. (10) — in terms of Ro, 
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The last equation, which is equivalent to i? = const = Ro, i.e. 

^[Ro,Q,P^]=0, (21) 

is automatically verified because is a quadratic function of Pr and hence dH/dPR is zero for circular orbits. 
We consider now a perturbation of the circular orbit defined by 

Pr = 5Pr , (22a) 

P^ = P°+ 5P9 , (22b) 

R = Ro+5R, (22c) 

w = Wo + 5a; . (22d) 

It is easy to verify that the Hamiltonian equations (17), when worked out at the linearized order, read as 

5Pr = -TTo 5R - po SPm, , (23a) 

<5P* = , (23b) 

6'R = aoSPR, (23c) 

6uj = po5R + To <5P* , (23d) 

where the coefficients, which depend on the unperturbed orbit, are given by 



dR"" 



[Po,0,P^], (24a) 

<Jo = ^^[Ro,Q,Pl] , (24c) 
oPr 

ro = 4^[Ro,0,P^] . (24d) 

By looking to solutions proportional to some e"^' one obtains some real frequencies, and therefore one finds stable 
circular orbits, if and only if 

Co = TTO (To > . (25) 

Using the Hamiltonian (7) we readily obtain 
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+ 



,.2 



, „ , rn- I 117 43 2 



-61 + 



-0' 

4777 



48 



325 2 

7r 

64 



31 O Q 



(26) 



This result does not look the same as our previous result (15), but this is simply due to the fact that it depends 
on the ADM radial separation Ro instead of the harmonic one ro- Fortunately all the material needed to connect 
Ro to ro with the 3PN accuracy is known [28] . In the case of circular orbits we readily find 



Ro =ro{l + 



„2 _4 



1 29 , 



.3 „6 



3163 21 
1680 ^ 



21 2 22 /ro\ 

32" -y^n^j 



^3 2 



(27) 



The difference between _Ro and ro is made out of 2PN and 3PN terms only. Inserting Eq. (27) into Eq. (26) and 
re-expanding to 3PN order we find that indeed our basic stability-criterion function Co comes out the same with 
our two methods. 

Finally let us give to the function Co an invariant meaning by expressing it with the help of the orbital frequency 
Wo of the circular orbit, or, more conveniently, of the frequency- related parameter 

^o=(^j • (28) 

From the inverse of Eq. (10) we readily obtain ro as a function of xo. This allows us to write the criterion for 

2 

stability as Co > 0, where Co = ^ 3 C0 admits the gauge-invariant form, which will be the same in all coordinate 

C 

systems. 



Co = 1 — 6 Xo -|- 14 z/ So + 



397 
2 



123 2 
16"" 



y - 14z/^ ) xl + 0{xt] 



(29) 



This form is more interesting than the coordinate-dependent expressions (15) or (26), not only because of its 
invariant form, but also because as we see the IPN term yields exactly the Schwarzschild result that the innermost 
stable circular orbit or ISCO of a test particle (i.e. in the limit f — > 0) is located at xisco = 1/6. Thus we find 
that, at the IPN order, but for any mass ratio v, 

IPN 1 



(30) 



One could have expected that some deviations of the order of u already occur at the IPN order, but it turns out 
that only from the 2PN order does one find the occurence of some non-Schwarzschildian corrections proportional 
to f. At the 2PN order we obtain 

3 



2PN _ 

a^isco — 



14i^ 



1- -i 1- 



(31) 



For equal masses this gives xflco — 0.187. Notice also that the eiTect of the finite mass corrections is to increase 

the frequency of the ISCO with respect to the Schwarzschild result (i.e. to make it more inward), and we find 
^ilco ~ i [1 + Ts'^ + ^i'^^)] ■ Finally, at the 3PN order and for equal masses u = j, we find that according to 
our criterion all the circular orbits are stable, and there is no ISCO. More generally, we find that at the 3PN order 
all orbits are stable when the mass ratio is > i^c where Vc — 0.183. 

Note that the above stability criterion Co gives an innermost stable circular orbit, when it exists, that is not 
necessarily the same as — and actually differs from — the innermost circular orbit or ICO, which is defined by the 
point at which the center-of-mass binding energy of the binary for circular orbits reaches its minimum value [41]. 
In this respect the present formalism, which is based on systematic post-Newtonian expansions (without using 
post-Newtonian resummation techniques like Pade approximants [42]), differs from some "Schwarzschild- like" 
methods such as the efi^ective-one-body approach [43] in which the ICO happens to be also an innermost stable 
circular orbit or ISCO. 

As a final comment, let us note that the use of a truncated post-Newtonian series such as Eq. (29) to determine 
the ISCO is a priori meaningful only if we are able to bound the neglected error terms. Furthermore, since we 
are dealing with a stability criterion, it is not completely clear that the higher-order post-Newtonian correction 
terms, even if they are numerically small, will not change qualitatively the response of the orbit to the dynamical 
perturbation. This is maybe a problem, and which cannot be answered rigorously with the present formalism. 
However, in the regime of the ISCO (when it exists), we have seen that xo is rather small, xo — 0.2 (this is 
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also approximately the value for the ICO computed in Ref. [41]), which indicates that the neglected terms in the 
truncated series (29) should not contribute very much, because they involve at least a factor Xq 0.002. On the 
other hand, we pointed out that in the limit z/ — >■ the criterion Co gives back the correct exact result, a^isco — |- 
This contrasts with the gauge- dependent power series (15) or (26) which give only some approximate results. 
Based on these observations, we feel that it is reasonable to expect that the gauge-invariant stability criterion 
defined by Eq. (29) is physically meaningful. 



Acknowledgements This article is based on the work [28] which has been done in collaboration with Bala R. 
Iyer. 
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